Network reconstruction from random phase-resetting 
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We propose a novel method of reconstructing the topology and interaction functions for a general 
oscillator network. An ensemble of initial phases and the corresponding instantaneous frequencies 
is constructed by repeating random phase-resets of the system dynamics. The desired details of 
network structure are then revealed by appropriately averaging over the ensemble. The method is 
applicable for a wide class of networks with arbitrary emergent dynamics, including full synchrony. 



o ■ 
(N : 



U 



> 
m 



X 



Complex networks of many interacting units found at 
all scales in nature are the subject of intense research in 
many scientific areas [l[ . Among the central issues in this 
field are the exploration and development of methods for 
determining the architecture of a network based on the 
observable data. Knowing the network structure helps in 
understanding its collective behavior, and indicates ways 
to engineer networks with the desired properties. For in- 
stance, it has been realized that inferring the topology 
of gene regulatory networks is crucial for completing our 
knowledge about the inner workings of cells [2|. Many 
real networks display modular and community structure 
that is essential for their functioning (3] and can be ex- 
tracted using a variety of methods [J] , as done for yeast 
metabolic network [5|. Reconstruction techniques often 
rely on examining the time-series of network dynamics 
that can reveal its interaction functions Q. The net- 
work topology can be detected by studying the inter- 
changes among its collective behaviors or investigating 
its response dynamics 0- Structural properties can be 
determined from various time-scales in the emergence of 
synchronization [§] , or by employing specific control the- 
ory methods 0. Recently proposed techniques involve 
noisy dynamical correlations between the nodes [loj , and 
even tackle models with non-equilibrium dynamics 11 1. 

However, existing reconstruction methods, that often 
use network models with single-node dynamics repre- 
sented by different types of oscillators [1^ , typically re- 
quire long time-series of dynamical data, or a certain level 
of complexity in the emergent dynamics 0, 01 • Since 
synchronization destroys the initial node-related infor- 
mation, detecting network topology in such cases is ex- 
tremely difficult. Some methods are applicable only to 
sparse or non-directed networks, often providing results 
with only a limited precision [lO|. 

In this Letter, we propose a novel method of re- 
constructing the topology and interactions of a general 
oscillator network. Our idea relies on repeatedly re- 
initializing the network dynamics (e.g. by performing 
random phase-resets), in order to produce an ensemble 
of the initial dynamical data. We design the quantities 
obtained by averaging this ensemble, whose values reveal 
the desired details of the network structure. Our method 
is applicable to any directed and weighted network, with 
general interaction functions and oscillator frequencies. 



and with arbitrary emergent dynamics, while avoiding 
the need for long time-series. 

In the context of phase-resets, one is typically inter- 
ested in the phase-resetting curves, which specify the sys- 
tem's response to weak external perturbations |l3i . They 
have been investigated both experimentally [ij] and the- 
oretically 15 , 3] , and been shown to contain properties 
relevant for determining network details such as cluster- 
ing [l7| . An algorithm for the estimation of neuron in- 
teraction and its stability based on phase-resets has been 
proposed [l3l. We here employ phase-resetting some- 
what differently, since our interest lies in the internal 
network interactions, rather than its response to stimuli. 
Contrary to [l3 |. we use phase-resets only as a natural 
way to re-initialize the dynamics of an oscillator network, 
without measuring the phase shifts occurring due to re- 
setting. 

Our model consists of N oscillators (nodes) , character- 
ized by their phases ipi £ [0, 27r) and natural frequencies 
uji. They are coupled pair- wise, via general 27r-periodic 
interaction functions fij with zero mean: 



(pi = OJi + 
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Models of this type include the famous Kuramoto model 
and its generalizations, widely used in theoretical stud- 
ies, as well as for describing specific experimental situa- 
tions [s, 12, 14]. The functions fij{(t>) are generally non- 
symmetric with respect to exchange of indices, and thus 
fully define the dynamical network (order of indices de- 
termines the direction of interaction). Network adjacency 
matrix given as Aij = sgn \fij \ specifies its topology. Dy- 
namics starts from a set of initial phases (i.p.) which we 
denote as = {ifi, . . . ifN)it = 0), chosen from a distri- 
bution p{ip) > normalized to (27r)^. The method is 
based on two assumptions: (i) wc arc able to arbitrarily 
re-initializc the network dynamics / times, by indepen- 
dently resetting the phases of all nodes to a new state ip; 
(ii) we arc able to measure all the values ipi, and all ini- 
tial instantaneous frequencies , each time the dynamics 
is re-initialized (for / = 1, . . . /). As we show in what fol- 
lows, the ensemble of data for / ^ 1 created under these 
assumptions yields the entire network structure. 

Introducing a 27r-periodic test-function g = g{^i — ^Pj) 
with zero mean, our aim is to compute the reconstruction 



2 



index 5*^., defined as: 
5.,[g] = (2^)-^ 

'[0,27r] 

Taking the functions fij in Eq.dl]) to be generally given 



dip giif; - ipj)ipj{if) . (2) 



by the Fourier series fij{4>) = '^ij smn(j)+blj cosn0, 
we obtain the following expression for Sij : 
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which is independent of the frequencies uji. The inte- 
gral over ipi vanishes unless i = k. This implies that 
if Aij = 0, the corresponding Sij — 0, independently of 
the choice of g. The non-zero entries of Sij directly reveal 
the presence of network links. In addition, matrix Sij de- 
tects the desired properties of the interaction functions 
for appropriately selected test-function g. In particular, 
using f/(0) = 26™*^ we obtain the Fourier harmonics of 
fij , which are the interaction parameters a|"^ and : 
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Computation of Sij for adequate g amounts for recon- 
struction of any dynamical network described by Eq. ([T]) . 
Depending on the properties of fij that are to be exam- 
ined, other choices of g are also possible. When deal- 
ing with the empirical interaction functions involving an 
unknown number of Fourier harmonics, a specifically de- 
signed g based on the experimental assumptions about 
fij might be useful. This result is largely independent 
of the frequencies Wi, the network's directedness, and the 
distribution p. In particular, it is also independent of the 
network's final dynamical state, whether dependent on 
p or not. However, a constant component in case of fij 
with non-zero mean cannot be detected, since its pres- 
ence is indistinguishable from the natural frequency lo. 

To practically implement our method, we need to con- 
vert the integral from Eq.(l2]) into an average involv- 
ing discrete non-uniformly distributed empirical data 
{fiYi^i {ipiYi^i- To that end, we represent the func- 
tion ipj(ip) using the kernel smoother Q{ip — ipi) as: 



(pj{ip) = 



The denominator is just the empirical density p{'p) = 
y^.i Q i'P — 'Pi) obtained via kernel distribution esti- 
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18j . Since the integration over ip already provides 
smoothing, we take Q{ip — tpi) — > 5{ip — t^j), and replace 
the Eq.® with a practical formula for Sij: 



Si] [g] = 



1=1 



(4) 



which is the average of empirical ipjg weighted by i. 

The most trivial way to obtain the ensemble {(Pilf^i 
would be to pick the values from a fixed distribution 
p(t^). Instead, we seek to mimic an experimentally fea- 
sible situation by performing / random phase-resets of 
the network dynamics, separated by the time interval 
T. Mathematically, this amounts to adding the term 
ELi Kusin{p, + au)S{t-lT) to the RHS of Eq.© 0. 
For each reset / and each oscillator i, we independently 
pick the kicking strength Kij from a zero mean Gaus- 
sian distribution with standard deviation K = 1, and 
the phase-shift ai^i uniformly from [0, 2tt). The ensemble 
is constructed by storing the phase values immediately 
after resets. The resulting artificially created ensemble 
has little in common with the natural distribution of 
phases, and can be considered as approximately inde- 
pendent. This is expressed by separability of p{ip) into a 
product of N one-dimensional distributions pi{ipi): 
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each of which wc determine from generated data using 
the kernel estimation method [31 . After each reset, the 
ensemble of <p is computed using a small time inter- 
val. The phase value prior to reset is of no importance, 
since our interest is not in the phase-resetting curves, 
but in modeling a realistic way to create the ensemble 
ip. The described procedure is quite similar to the recent 
experimental implementation of the randomized phase- 
resetting of epileptic neurons aimed at their transient 
desynchronization j20l |. In these experiments, however, 
the problem of simultaneous read-out of phases and fre- 
quencies remains a challenge. 

We now illustrate our theoretical findings through nu- 
merical simulations on simple network examples, com- 
puting the reconstruction index Sij as described above. 
Consider a simple network with TV = 4 oscillators shown 
in Fig[TJ We pick the natural frequencies at random from 




FIG. 1: 4-node network used for illustrating our method. 



uji e [—1,1]. The interaction functions fij are defined 

for linked node pairs by randomly choosing Q-ij^ ^Ij' & 

[—1,1], while taking a^-J'' = b\^j^ = for n > 2. Since such 
a network typically does not synchronize, our approxi- 
mation of independent i.p. after resetting is appropriate. 
We take g = 2e^'^ and compute Sij from an ensemble of 
/ = 10"* i.p. to obtain the numerical approximations of 



-,(1) 



and b[^'^ via Eq.(j4]). In Figl2]we compare the numeri- 
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cal a-j'' and b[^^ (crosses) with the actual values (circles). 
All values display a very good agreement for both linked 
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FIG. 2: Reconstruction of the network from Fig[l] Circles: 
actual parameter values, crosses: numerically obtained values 
for 7 = 10''. Left: 
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right: fo'J', for each node pair i 



node pairs (different from zero) and non-linked node pairs 
(zero). We have not only revealed the adjacency matrix 
Aij, but also found the interaction parameters a^j'' and 

bl^\ thus reconstructing the entire dynamical network. 

Below we discuss the limitations of our method. If the 
available data ensemble / is too small, the statistics is 
poor and the obtained network characteristics have large 
uncertainties, which typically decrease as ^ To 
illustrate this, in Fig|3^ we present the numerical values 
of parameter a[j'', computed for network in Fig |1] using 
the ensemble of i.p. of size /. While the distinction 
between links and non-links can already be seen for / ^ 
10'^, for good approximation one needs / > 10* (as done 
in Figl2]). For higher Fourier harmonics, the convergence 
is gradually slower, but maintains the same properties. 

Another limitation is related to the validity of our in- 
dependence assumption for the ensemble of i.p. which is 
expressed by the separability of distribution p{ip) Eq.dS]). 
This heavily depends on the network's dynamical regime 
and the resetting strength. For a full synchrony and weak 
kicking, the reset state is expected to be strongly corre- 
lated, whereas for chaotic dynamics and strong resets, the 
independence assumption is essentially correct. To study 
this, we consider again the network from Figll] but now 
we fix all frequencies to Wi = 1, and take all interactions 



to be attractive a. 
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(Kuramoto- type model 
with identical oscillators). We apply random kicking as 
described above after allowing the network to synchro- 
nize (t 3> t synch), but this time with a variable standard 
deviation of kicking strength < K < 10. For each value 
of K we create an ensemble of / = 10* i.p., and use it to 
compute ajj' as done previously. In Fig[3lD we show the 

reconstructed values of a,-^^ for links and non-links in re- 
lation to K. Sufficiently strong kicking {K > 5) succeeds 
in destroying the network's synchrony and generating the 
independent i.p., from which a good approximation of 
a-j^ is computed. Moderate kicking ~ 1 applied pre- 
viously are now insufficient. This furthermore depends 
on the relation between r and t synch- if i$ t synch (fre- 
quent resets) the separability of p is easier to achieve. 



Too strong kicking can also induce correlations in (p, re- 
gardless of dynamical regime and r. However, note that 
p can be estimated using the techniques more elaborate 
than simple one-dimensional kernels [l8| . which can in 
principle yield a good estimate even in the non-separable 
case. On the other hand, phase-resetting is potentially 
not the only mechanism of obtaining the ensemble ip\ re- 
call that our theory with a known p{ip) works equally well 
for any case, including full synchrony and inseparability. 
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FIG. 3: (color online). Numerical values of a'j' for network 
in Fig[l] for links (cyan/gray) and non- links (black), (a) com- 
puted from ensemble of / i.p. (cf. Fig[2]). (b) computed from 
/ = 10* for network with attractive interactions, and with re- 
setting done at synchronous state using kicking strength K. 
(c) computed from / = 10^ for network with attractive inter- 
actions where only spikes (</? = 0) are observable, in relation 
to coupling strength e (see text for details). 



Adding noise terms to RHS of Eq.([T|) does not formally 
change the derivation of our main result, rendering our 
theory valid in the presence of noise. However, in the 
light of discussion above, noise will have an effect on the 
performance of method: additional uncertainty due to 
larger fluctuations of the estimated (p require larger en- 
sembles to achieve the desired precision. On the other 
hand, noise may play a constructive role by destroying 
the undesired correlations within ip, and thus facilitating 
the separability of p. 

While the experimental techniques for measuring ip 
arc already in use [l3 |. in a potential realistic applica- 



4 



tion of our method a problem may arise in relation to 
the measurement of (p. The entire cycle of a real os- 
cillator is often not accessible; instead, one can observe 
only a single event per period (e.g. a spike produced by 
a neuron). In such cases, one is forced to estimate the 
instantaneous frequencies relying solely on the time in- 
tervals between the spikes. To illustrate this, we consider 
again the system studied in FiglSlD, but now we replace 
a-j^ with ea-j^. The parameter e (coupling strength) con- 
trols the ratio between the oscillation time-scale (period) 
and the interaction time-scale (synchronization). Rather 
than computing instantaneous (p after each reset, we ob- 
serve only the event of an oscillator passing through the 
phase value = (spike), and estimate both and (p 
from the first two spikes observed after resetting. We 
then reconstruct the values of ai^ using the ensemble of 
/ = 10^ i.p. as done before (strong resetting is applied 
immediately after the spikes are recorded). The results 
shown in FiglSt have a clear physical interpretation: for 
too small coupling e < 0.03 the links can not be revealed 
since the interaction is too weak. For too large e > 0.4 
the two time-scales are too close, and the detection is 
again impossible since the distribution of phases changes 
significantly over a period. However, between these ex- 
tremes, there is a range of coupling around e 0.1 where 
the two time-scales are well separated allowing a reliable 
reconstruction. This shows that with an adequately big 
ensemble our method works even if the entire oscillator 
cycle is not accessible: errors in the estimation of ip and 
If play a role similar to the noise. The method fails in 
the case of too strong coupling, similarly to the case of 
too weak resetting after synchronization (cf. FigI3]3). 

In conclusion, we proposed a method of reconstruct- 
ing oscillator networks by repeating random phase-resets, 
applicable to a general network irrespectively of the dy- 
namical regimes (the feasibility of such resetting has been 
recently demonstrated for neural tissue [l^l). Our the- 
ory emphasizes the importance of the transient dynam- 
ics in the context of network reconstruction, thus com- 
plementing the available techniques that rely on time- 
series recorded in final stationary state. Our theoreti- 
cal model can be straightforwardly generalized to other 
models beyond Eq.dT]). If the couplings depend on two 
phases in a more general way, or depend on more than 
two phases, one should use more elaborate test-functions 
(e.g. in a form of general complex exponentials); how- 
ever, even a theoretical description of such networks is 
already a challenge. For high-dimensional oscillators only 
a single scalar might be observable: our method can 
still be applied through the appropriate transformation 
to phases [y] . Another generalization regards the recon- 
struction of sub-networks, in the case that only infor- 
mation on some nodes is accessible. The problem here 
is to infer the distribution of i.p. for the non-accessible 
nodes. Finally, a real experimental situation may involve 
a network whose dynamics cannot be reset for all nodes 



simultaneously, which renders the independence assump- 
tion invalid. This is a much more challenging, although 
very realistic case that requires additional study. 
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